# All the code necessary to make those cool maps of voting by state, income, and ethnicity

library (arm)
source("FUN.R")


naes04 <- read.dta ("naes04.dta")###
year=2004

dataset <- naes04

complete <- !is.na(dataset$stnum) & !is.na(dataset$ethn) & !is.na(dataset$inc) & !is.na(dataset$favvoucher)

# (1) Some info that we'll need

n.state <- 51
n.inc <- 5
n.eth <- 4
region <- c(3,4,4,3,4,4,1,1,5,3,3,4,4,2,2,2,2,3,3,1,1,1,2,2,3,2,4,2,4,1,1,4,1,3,2,2,3,4,1,1,3,2,3,3,4,1,3,4,1,2,4)
state.abb.long <- c (state.abb[1:8],"DC",state.abb[9:50])
state.name.long <- c (state.name[1:8],"Washington, D.C.",state.name[9:50])

wmean <- function (a, w=rep(1,length(a)), subset=rep(TRUE,length(a))){
  keep <- !is.na(a) & !is.na(w) & !is.na(subset) & subset
  sum ((w*a)[keep])/sum(w[keep])
}

votes <- read.dta ("state vote and income, 68-00.dta")
income2000 <- votes[votes[,"st_year"]==2000, "st_income"]
state.income <- c (income2000[1:8],NA,income2000[9:50])
z.state.income <- rescale (state.income)

vote2004 <- read.table ("2004data.dat", header=TRUE)
vote2004.long <- rbind (vote2004[1:8,], c(184,19), vote2004[9:50,])
rvote04 <- 1 - as.vector (vote2004.long[,1]/rowSums(vote2004.long))

votes08 <- read.csv ("2008ElectionResult.csv")
obama08 <- votes08[,"vote_Obama"]
mccain08 <- votes08[,"vote_McCain"]
rvote08 <- mccain08/(obama08+mccain08)

namelist <- list (state.abb.long,
                  c("0-20K","20-40K","40-75K","75-150K",">150K"),
                  c("White Catholics","White Evangelical Protestants","White Non-Evang. Protestants","White Other or No Religion", "Blacks", "Hispanics", "Other Races"),
                  c("unweighted","weighted"))
namelist2 <- list (namelist[[1]], namelist[[2]])
namelist3 <- list (namelist[[1]], namelist[[2]], namelist[[3]])

# (2) Load in the Annenberg data

y <- dataset[complete,"favvoucher"] ###
y <- ifelse (y==1|y==2, 0, ifelse (y==4|y==5, 1, NA)) ###
inc <- dataset[complete,"inc"] ###
eth <- dataset[complete,"ethn"] ###
st <- dataset[complete,"st"] ###
state <- rep (NA, length(y))
for (i in 1:51){
  state[st==state.abb.long[i]] <- i
}

relig <- dataset[complete,"relig"] ###
relig[is.na(relig)] <- 7
bagain <- dataset[complete,"bagain"] ###
bagain[is.na(bagain)] <- 0
releth <- ifelse (relig==2 & eth==1, 1,
            ifelse (relig==1 & eth==1 & bagain==1, 2,
              ifelse (relig==1 & eth==1 & bagain==0, 3,
                ifelse (eth==1, 4,
                  ifelse (eth==2, 5,
                    ifelse (eth==3, 6, 7))))))
n.releth <- max (releth)

pop.weight <- rep (1, length(y))

# Make the 3-way dataset
ybar.weighted <- array (NA, c(n.state,n.inc,n.releth), dimnames=namelist3)
n <- array (NA, c(n.state,n.inc,n.releth), dimnames=namelist3)
design.effect.by.cell <- array (NA, c(n.state,n.inc,n.releth), dimnames=namelist3)
for (i in (1:n.state)){
  if (!(state.abb.long[i] %in% c("AK","DC","HI"))){
    for (j in 1:n.inc){
      for (k in 1:n.releth){
        ok <- state==i & inc==j & releth==k
        ybar.weighted[i,j,k] <- wmean (y, pop.weight, ok)
        n[i,j,k] <- sum (ok)
        design.effect.by.cell[i,j,k] <- 1 + var(pop.weight[ok]/mean(pop.weight[ok]))
      }
    }
  }
}
ybar.weighted[is.nan(ybar.weighted)] <- NA
design.effect <- wmean (design.effect.by.cell, n, n>1)
n.eff <- n/design.effect

# (3) Read in the 3-way table of voter turnout and adjust it for religion

load ("prop.voters.Rdata")
prop.voters0 <- prop.voters
prop.voters <- NA + ybar.weighted
prop.voters[,,5:7] <- prop.voters0[,,2:4]
for (i in 1:n.state){
  for (j in 1:n.inc){
    ok <- state==i & inc==j
    counts <- rep (NA, 4)
    for (k in 1:4){
      counts[k] <- sum (ok & releth==k)
    }
    prop.voters[i,j,1:4] <- prop.voters0[i,j,"Whites"]*counts/sum(counts)
  }
}

# (4) Fit an lmer model on the 3-way dataset

# Function to create an index variable by combining two or more separate
# index variables, each of which must take on positive integer values
cross <- function (...){
  a <- list (...)
  n.factors <- length (a)
  multiplier <- 1
  b <- 0
  for (j in 1:n.factors){
    b <- b + multiplier*a[[j]]
    multiplier <- multiplier * round (10^(ceiling(log10(max(unique(a[[j]]))))))
  }
  return (b)
}

# Put the data vectors together
count <- 0
state.v <- rep (NA, n.state*n.inc*n.releth)
inc.v <- rep (NA, n.state*n.inc*n.releth)
releth.v <- rep (NA, n.state*n.inc*n.releth)
for (k in 1:n.releth){
  for (j in 1:n.inc){
    for (i in 1:n.state){
      count <- count + 1
      state.v[count] <- i
      inc.v[count] <- j
      releth.v[count] <- k
    }
  }
}
z.inc.v <- rescale (inc.v)
z.state.inc.v <- z.state.income[state.v]
rvote04.v <- rvote04[state.v]
region.v <- region[state.v]
inc.region.v <- cross (inc.v, region.v)
inc.state.v <- cross (inc.v, state.v)
inc.releth.v <- cross (inc.v, releth.v)
releth.region.v <- cross (releth.v, region.v)
releth.state.v <- cross (releth.v, state.v)

n.eff.v <- as.vector (n.eff)
ybarw.v <- as.vector (ybar.weighted)
ybarw.v[n.eff.v==0] <- .5
data.v <- cbind (ybarw.v*n.eff.v, (1-ybarw.v)*n.eff.v)

threeway.fit <- glmer (data.v ~ z.inc.v*z.state.inc.v + z.inc.v*rvote04.v + (1 + z.inc.v | region.v) + (1 + z.inc.v | state.v) + (1 | inc.region.v) + (1 + z.inc.v| releth.v) + (1 | inc.releth.v) + (1 | inc.state.v) + (1 | releth.region.v) + (1 | releth.state.v) + (1 | inc.v) , family=quasibinomial(link="logit"))
display (threeway.fit)

fit <- rep (NA, n.state*n.inc*n.releth)
fit[!is.na(ybarw.v)] <- fitted (threeway.fit)

theta.hat.v <- fit

# Put the estimated cell means back as a three-way array
theta.hat.unshifted <- array (theta.hat.v, c(n.state,n.inc,n.releth),
                              dimnames=namelist3)

# (5) For estimating opinion on a survey question, no renormalization needed

theta.hat <- theta.hat.unshifted

# (6) Print the estimates and make the maps

# Print everything
pfround (prop.voters, 2)
pfround (ybar.weighted, 2)
pfround (theta.hat, 2)


#  year <- coolmaps.info$year
#  theta.hat <- coolmaps.info$theta.hat
#  prop.voters <- coolmaps.info$prop.voters
#  ybar.weighted <- coolmaps.info$ybar.weighted

# Create the totals, summing over all ethnic groups
ybar.weighted.2way <- array (NA, c(n.state,n.inc), dimnames=namelist2)
theta.hat.2way <- array (NA, c(n.state,n.inc), dimnames=namelist2)
n.eff.2way <- array (NA, c(n.state,n.inc), dimnames=namelist2)
for (i in 1:n.state){
  for (j in 1:n.inc){
    ybar.weighted.2way[i,j] <- wmean (ybar.weighted[i,j,], prop.voters[i,j,])
    theta.hat.2way[i,j] <- wmean (theta.hat[i,j,], prop.voters[i,j,])
    n.eff.2way[i,j] <- sum (n.eff[i,j,])
  }
}

# Make the maps for all voters and whites only
inclabel0 <- c ("under $20,000",
                "between $20,000 and $40,000",
                "between $40,000 and $75,000",
                "between $75,000 and $150,000",
                "over $150,000")
inclabel1 <- c("Income under $20,000", "$20-40,000", "$40-75,000", "$75-150,000", "Over $150,000" )

# Function for a blank plot
blankplot <- function (words, cex=1){
  plot (c(0,1), c(0,1), xlab="", ylab="", xaxt="n", yaxt="n", bty="n",type="n")
  text (.5, .5, words, cex=cex)
}

prop.voters.sum <- array (NA, c(n.state, n.releth))
for (i in 1:n.state){
  for (k in 1:n.releth){
    prop.voters.sum[i,k] <- sum (prop.voters[i,,k])
  }
}

# Make the map  
png (paste ("vouchermaps", year, ".png", sep=""), height=1000, width=1000)
par (mar=c(0,0,2,0), mfcol=c(8,6), tck=-.01, mgp=c(1.5,.5,0), oma=c(7,0,5,0))
national.avg <- wmean(theta.hat,prop.voters)
for (k in 0:7){
  blankplot (c("All voters", "White\nCatholics", "White evangelical\nProtestants", "White non-evang.\nProtestants", "White other/\nno religion", "Blacks", "Hispanics", "Other races")[k+1], cex=2)
}
for (j in 1:n.inc){
  for (k in 0:7){
    if (k==0){
      statemaps (orangegreen2 (-(theta.hat.2way[,j]-national.avg), lo=-.25, hi=.25))
      mtext (inclabel1[j], 3, cex=1.2, line=1)
    }
    else {
      statemaps (ifelse (prop.voters[,j,k]>.01, orangegreen2 (-(theta.hat[,j,k] - national.avg), lo=-.25, hi=.25), "white"))
    }
  }
}
mtext (paste (year, ":  State-level support (orange) or opposition (green) on school vouchers, relative to the national average of ", round (100*national.avg), "% support", sep=""), side=3, outer=TRUE, line=3, cex=1.3)
mtext (paste ("Orange and green colors correspond to states where support for vouchers was greater or less than the national average.\nWhere a category represents less than 1% of the voters of a state, the state is left blank.", sep=""), side=1, outer=TRUE, line=5, cex=1)
dev.off()

# (7) Make graphs for the 50 states

for (k in 1:7){
  png (paste ("voucherplots", year, " ", namelist[[3]][k], ".png", sep=""), height=1000,width=1000)
  par (mar=c(2,2,1,0), tck=0, mgp=c(1.5,.5,0), oma=c(0,0,4,0))
  graph.dims <- c(7,7)
  par (mfrow=graph.dims)
  sort.state <- rev (order (rvote08))
  count <- 0
  for (i in sort.state){
    if (!(state.abb.long[i] %in% c("AK","HI","DC"))){
      count <- count + 1
      if (prop.voters.sum[i,k] > .1){
        plot (c(1,5), c(0,1), xlab="", ylab="", xaxt="n", yaxt="n", type="n",
              yaxs="i")
        if (count%%graph.dims[2]==1)
          axis (2, c(.02,.5,.95), c("0","50%","100%"), cex.axis=1.1)
        if (count > (48-graph.dims[2]))
          axis (1, c(1.2,3,4.8), c("poor","mid","rich"), cex.axis=1.2)
        abline (.5, 0, col="gray", lwd=.5)
        fit <- theta.hat[i,,k]
        pts <- ybar.weighted[i,,k]
        neff <- n.eff[i,,k]
        se <- sqrt (fit*(1-fit)/neff)
        text (3, .9, state.name.long[i], cex=1.3)
        for (j in 1:n.inc){
          lines (rep (j, 2), pts[j] + c(-1,1)*se[j], lwd=.5)
          points (j, min (.98, max (.02, pts[j])), pch=20, cex=1.5)
        }
        lines (1:n.inc, fit, lwd=2)
      }
      else {
        plot (c(1,5), c(0,1), xlab="", ylab="", xaxt="n", yaxt="n", bty="n", yaxs="i", type="n")
        text (3, .9, state.name.long[i], cex=1.3)
      }
    }
  }
  mtext (paste (year, ":  Raw and estimated support for school vouchers within each state among", namelist[[3]][k]), outer=TRUE, cex=1, side=3, line=1)
  dev.off()
}



coolmaps.make()
